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A dynamic localization model 
with stochastic backscatter 

By D. Carati &: S. Ghosal 

1. Motivation and objectives 

1.1 The dynamic localization ’procedure 

The modeling of subgrid scales in large-eddy simulation (LES) has been rational- 
ized by the introduction of the dynamic localization procedure (Ghosal et al. 1993, 
1994). This method allows one to compute rather than prescribe the unknown co- 
efficients in the subgrid-scale model. Formally, the LES equations are supposed to 
be obtained by applying to the Navier-Stokes equations a “grid filter” operation 
defined as: 


«/>(*) = j y dy G (^) V’(y). (1) 

where G is a kernel damping the fluctuations with a characteristic length shorter 
than A. The resulting equations (here we only consider incompressible flows) 

dt u i H* dj(ujUi) = i/oV 2 iij — dip Oj T t j , (2) 

contain an unknown “subgrid stress” tensor t 1} that needs to be modeled: 

Tij - u~Uj - Ui Uj. ( 3 ) 

Though the subgrid stress itself is unknown, an identity between subgrid stresses 
generated by different filters has been derived (Germano et al. 1991) and is the 
basic ingredient of the dynamic procedure: 

Lij = Tij — T{j , ( 4 ) 

where Lij = uiuj —Ui Uj is the Leonard tensor and Tij = u, Uj — u, Uj is the subgrid 
stress tensor generated by a second filter defined by: 

^(x) = f dy G{ - ^ — ) V(y)- ( 5 ) 

Jv A 

Here G is a kernel damping the fluctuations with a characteristic length shorter than 
A. It will be referred to as the a test filter”. If models are used for these 

quantities, the difference between the right- and the left-hand sides of relation (4), 

Eij = Lij + T% -T% ± 0, 


(6) 
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may be used as “quality indicators” for the subgrid-scale models. In practice, if the 
models contain a small number of unknown parameters (like in the Smagorinsky 
(1963) model t, } — ^SijTkk = - 2CA 2 |,S|Sjj), the dynamic procedure proposes to 
determine these parameters by minimizing the quantity: 

F\C] = J v dy Eij[y; C] E tJ [y; C). (7) 

This procedure partially removes the arbitrariness inherent to modeling in LES. 
However, the success of the dynamic procedure still strongly depends on the quality 
of the model for which it is implemented. 

1.2 Dynamic localization model with k-equation 

This model is motivated by the following considerations. When no constraint 
is imposed on the Smagorinsky coefficient C, the minimum of the functional E[C) 
is achieved for a field C, which can be locally negative. In this case, the model 
exhibits local reverse energy transfer. This is one of the simplest adaptations of 
the standard Smagorinsky expression to allow for a variable C and backscatter 
(Piomelli et al. 1991). However, it leads to some difficulties. A negative eddy 
viscosity generates an exponential amplification of local disturbances instead of the 
traditional exponential damping. The resulting backscatter is an “auto-catalytic” 
phenomenon which does not correspond to the real physics of reverse energy transfer 
in turbulent flows. As a consequence, unphysical instabilities in the LES equations 
have been observed when the coefficient C in the Smagorinsky model is determined 
by an unconstrained (no positivity required) variational procedure. The backscatter 
appears to be unsaturated and the model is unusable. However, Ghosal et al. (1993, 
1994) have stressed that the reverse flow of energy must be quenched at the latest 
when all the subgrid-scale energy has been removed, and they have proposed to use 
the alternate representation 


i/, = 2C'A k 1 ' 2 (8) 

instead of the Smagorinsky scaling. Here, k represents the subgrid-scale energy for 
which a separate transport equation is needed. The basic DLM(k) equations axe 
given in (Ghosal et al. 1993, 1994). It can be shown that this model is stable. This 
approach involving the subgrid-scale kinetic energy is the first self-consistent model 
in which backscatter is accounted for in the framework of the dynamic procedure. 

The combination of a locally negative transport coefficient and a saturation pro- 
cess is reminiscent of some instabilities in complex fluids. For example, phase sep- 
aration in multicomponent mixtures can be described by instabilities created by 
a negative diffusion coefficient and saturated by surface tension effects (which are 
usually modeled by “hyperdiffusivity” terms). Roughly speaking, the DLM(k) pic- 
ture for backscatter is similar. At some locations in the fluid, the eddy viscosity 
becomes negative. In a first stage, this generates an instability characterized by an 
exponential amplification of the local disturbances. In a second stage, a saturation 
process arrests the further growth of the instability. Later, the rapid changes in 
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the turbulent velocity are likely to modify the local conditions, and the viscosity 
is expected to go back to positive values. This process does not explicitly take 
into account the possible stochastic nature of backscatter, but is not incompatible 
with a “molecular representation” of the small eddies. Indeed, the large diversity 
of small-scale eddies suggests that the turbulent fluid should behave more like a 
very complex (in a rheological sense) medium, and one should not expect the eddy 
viscosity to remain positive at every space time point. 

Although preliminary tests of this model have been satisfactory, the use of a neg- 
ative eddy viscosity to describe backscatter is probably a crude representation of 
the physics of reverse transfer of energy. Indeed, the model is fully deterministic. 
Knowing the filtered velocity field and the subgrid-scale energy, the subgrid stress 
is automatically determined. Obviously, this is only an approximation. It is very 
unlikely that the small scales influence the large scale evolution only through k. 
This is nevertheless an improvement when compared to the traditional Smagorin- 
sky model in which no information from the small scales is included. However, we 
know that the LES equations cannot be fully deterministic since the small scales are 
not resolved. This stems from an important distinction between equilibrium hydro- 
dynamics and turbulence. In equilibrium hydrodynamics, the molecular motions 
are also not resolved. However, there is a clear separation of scale between these 
unresolved motions and the relevant hydrodynamic scales. The result of molecular 
motions can then be separated into an average effect (the molecular viscosity) and 
some fluctuations. Due to the large number of molecules present in a box with 
size of the order of the hydrodynamic scale, the ratio between fluctuations and the 
average effect should be very small (as a result of the “law of large numbers”). For 
that reason, the hydrodynamic balance equations are usually purely deterministic. 
In turbulence however, there is no clear separation of scale between small and large 
eddies. In that case, the fluctuations around a deterministic eddy viscosity term 
could be significant. An eddy noise would then appear through a stochastic term 
in the subgrid- scale model and could be the source of backscatter. Some existing 
models have already represented reverse energy transfers by random terms. For 
example, a random eddy force derived from the eddy damped quasi-normal Marko- 
vian approximation has been used with some success in LES of isotropic turbulence 
by Chasnov (1991). This idea has been extended to boundary layers by Mason 
& Thomson (1992) and a similar approach has also been used by Leith (1990) to 
study LES of mixing layers. However, all these stochastic models contain an arbi- 
trary parameter that must be tuned to obtain satisfactory results. Here we present 
an alternative subgrid-scale model in which the dynamic procedure is combined 
with a stochastic representation of backscatter. Following the dynamic procedure, 
no arbitrary parameter will be introduced in the model. Such a model represents a 
more traditional picture (Kraichnan, 1976; Leslie & Quarini, 1979) of backscatter 
than a negative eddy viscosity based model. However, it must be stressed that 
the true energy transfers between small and large scales are probably much more 
complex than that described by either an eddy viscosity or an eddy noise formal- 
ism. Both these models probably remain rather crude approximations to the real 
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physical process. 

2. Accomplishments 

2.1 Stochastic dynamic localization model 

The DLM(k) proposed by Ghosal et ai (1993, 1994) accounts for backscatter 
through a purely deterministic eddy viscosity. Let us now adopt a different point of 
view and assume that backscatter may be represented by a stochastic forcing term. 
At grid level, the proposed model is: 

d } T tJ = d } (C M + /., (9) 

where ftij = —2 A 2 \S\ Sij corresponds to the standard Smagorinsky model and fi 
is an eddy force. For the sake of simplicity we will choose the simplest temporal 
behavior for f by assuming that the eddy force is a white noise process. The general 
form of its two-point, two-time correlation is then given by: 

(fi(r>t)fj(r',t')) = A 2 (r,t)Hij(r — r')6(t — t'), (10) 

where (r, t) and (r', t f ) are two space-time points. The operator (...) will denote the 
averaging over all possible realizations of the random force conditioned on a given 
velocity field u(r, t ). The functions Hij characterizing the forcing correlation will be 
discussed later. We only assume that the prefactor A 2 is chosen so that i?j,(0) = 1. 
In what follows, /, is supposed to be divergence free (this can always be ensured 
by suitably modifying the pressure term). The choice of a solenoidal force is not 
essential but it simplifies the following discussion because the pressure then does 
not involve the random fields used to model the stochastic force. In some sense, the 
white noise process can be seen as the “most stochastic” choice. Thus, comparison 
between the stochastic dynamic localization model, DLM(S) defined by (9), and the 
DLM(k) should show what are the respective advantages (if any) of stochastic and 
deterministic models for backscatter. 

It should be noted that the DLM(S) only models the divergence of the subgrid- 
scale stress. This is justified because the divergence is the only quantity needed in 
the LES equations. Also, the introduction of a stochastic force is much easier in the 
formulation (9). Thus, the quantity djEij should be used in the dynamic procedure 
instead of Eij . However, this would lead to major difficulties: The unknown quan- 
tities A and C would be determined by stochastic, integro-differential equations. 
The resolution of such equations would dramatically reduce the performances of 
the model. To avoid these problems, we propose to base the minimization proce- 
dure on the quantity (Eij) instead of djEij where the average is performed over all 
the possible realizations of the random noise /, conditioned on a fixed velocity field 
U (r,0- This is a convenient approximation which results in the following simplifi- 
cations. First, the error tensor (Eij) is deterministic and totally independent of the 
random forces: 


(E{j) — Lij + C fiij — C 


(ii) 
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where aij = —2 A 2 |5| 5,j. Also, C is now determined by minimizing J y dy(Eij)(Eij) 
which is exactly the same variational problem as in the deterministic models. Fi- 
nally, since the DLM(S) is supposed to model backscatter by the eddy force, it is 
natural to assume that the Smagorinsky term is purely dissipative. The parame- 
ter C has then to be determined by the constrained dynamic localization model, 
DLM(+) (see Ghosal et al. 1993, 1994). 

Since the force disappears from this first part of the dynamic procedure due 
to the conditional averaging over all realizations of the random force, we need an 
extra-relation for determining its amplitude. It can be obtained by noting that, 
even though the average effect of the force vanishes in the equations for the mean 
momentum, it will lead to a finite effect in the energy balance equation. Two 

equivalent energy equations may be obtained for the quantity E = UiUi/ 2: 

DtE — . . . — u^dj (C7 Qfjjf -j- P (12a) 

dtE — . . . — U{dj [c 0ij + Lij + p (126) 

where . . . stands for the viscous and inertial terms that axe identical in both these 
equations. The pressure terms P and p are determined to keep the velocity diver- 
gence free at grid and test level. The quantities £p and £j represent the energy 

input in the system respectively by Fi (the test level eddy force) and /, (the filtered 
grid level eddy force). The difference between the right-hand sides of Eqs. (12a) 
and (12b) 


Z - £ F - £j - g ^ 0 (13) 

plays exactly the same role for the energy transfer as the quantities Eij for the 
subgrid-scale stress. Here g is a known quantity ( C has been determined by the 
DLM(+)) given by 

9 = Uid } (c a, } + P 6ij - C fcj - Lij - p Sij'j . (14) 

The minimization of the quantity Z — f y dr (Z) 2 can now be used as a variational 
determination of the parameters that enter the model for the eddy force. At this 
point, little has been said about the statistical characteristics of the stochastic force 
itself. The variational problem presented here could be used together with a wide 
variety of choices for the eddy force. 

Let us now discuss some additional assumptions that will greatly simplify the 
DLM(S) equations. From a computational point of view, it will be very convenient 
to consider the limiting case for which /, at different grid points may be assumed 
to be completely uncorrelated. This avoids the non-trivial problem of generating 
random numbers with complex spatial correlations. This is also physically plausible 
since the eddy force is assumed to model random phenomena due to structures 
smaller than the mesh grid. Thus, the function H t j will be assumed to be negligible 
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for distances |r — r'| larger than the mesh grid. Here we also assume that the 
probability distribution function of /, corresponds to a Gaussian process. In that 
case, a force characterized by the correlations (10) leads to an average energy input 

£ f = ±A 2 (r,t)H ii (0) = ±A 2 (r,t). (15) 

Let us now show how the dynamic procedure can be used to determine this energy 
input. The variational formulation presented above does not directly involve £/ . 
However, we can easily relate it to £/ using Kolmogorov type ideas about energy 
transfer in the inertial range. Indeed, all the models used in LES - and model (9) 
is not an exception - are based on such arguments. Thus, the transfer rate is 
supposed to be independent of the filter width (in the inertial range), and we may 
assume that the backscatter rate has the same property (£/) = (£f)- Moreover, the 
test-filter acts like a local averaging of the random numbers /’s which are almost 
uncorrelated (the two-point correlation is assumed to decay rapidly) and which 
vanish on the average. Thus, / is much smaller than /, and (£y) can be neglected 
when compared to (£/). It can then be assumed that (£^) (£/) = (£f)> Thus, 

relation (13) gives 


(Z) = (£ f )-g. (16) 

Let us now consider a simple model for the stochastic force: 

f, = Pij(A e,) (17) 

where A is a (dimensional) coefficient which plays a role similar to the Smagorinsky 
coefficient C. The operator P u = tfy - V~ 2 V,Vj takes out the divergence of the 
vector At y Here e; are random numbers for which the probability distribution 
function is supposed to be Gaussian and isotropic: 

( e i( r , <)) = 0, (18a) 

<e,-(M)e>(r',f')) = 1 6 tj S(t - t') 6 r , r ., (186) 

Here, we focus on the discrete equations and consequently we have used the Kro- 
neker symbol. In agreement with the arguments leading to relation (16), the two- 
point correlation vanishes for r ^ r'. Comparison of (17) and (18) with (10) shows 
that Hij( r,r') = PikPjkt r,r'/ 2 and A 2 = 2A 2 /3. Thus, the energy input is: 

(£ f ) = ^A 2 = ^A 2 . (19) 

The variational problem can then be used to determine A by minimizing 

T~ S ) ■ 


( 20 ) 
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FIGURE 1. Time evolution of spectra in decaying isotropic turbulence : 

DLM(k); +: DLM(+); : DM; : DLM(S); : no model; •: experi- 

ment (t=1.55); ■ : experiment (t=2.70). 


Following the method used by Ghosal et a l. (1993, 1994) this leads to 


A 2 = 3 [g] + , (21) 

where [x]+ = (x + |x|)/2. Thus, if we consider that the random numbers e, are 
uncorrelated for different grid points, the DLM(S) is defined by (9) and (17) in 
which C has to be determined by the DLM(-|-) and the forcing amplitude A is 
given by the explicit relation (21). 

2.2 Results 

The model described in the previous section has been extensively tested and com- 
pared to other models (the original dynamic model: DM; the constrained dynamic 
localization model: DLM(+), and the DLM(k)) for forced and decaying isotropic 
flows. We will not discuss in detail the conditions of these simulations which are the 
same as those presented in previous reports (Ghosal et al. 1993 and Ghosal 1993). 

Fig. 1 shows the result of a simulation performed using 48 grid points in each 
direction. This seems to be the smallest simulation that would be consistent with 
the condition implicit in the idea of LES, viz., that the subgrid scales should carry 
significantly less energy than the resolved scales. 
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Figure 2. Prediction of Kolmogorov’s 5/3 law and the Kolmogorov constant in 

forced isotropic turbulence at steady state : DLM(k) - +• DLMf+V • 

DM; : DLM(S). 

All four models predict decays in the resolved energy that are in good agreement 
with the experiment (Comte-Bellot and Corrsin 1971). 

In the asymptotic self-similar regime the energy decays as a power law E ~ t~ a . It 
is not clear that such a self similar regime is reached in the present experiment since 
only three experimental points are available. However, the three experimental points 
almost lie on a straight line on a log-log plot. The decay exponent is thus estimated 
to be a r; 1.20. A least-square fit to the LES data yields the values a — 1.27, DM - 
a = 1.21, DLM(+); a = 1.28, DLM(k); a = 1.17, DLM(S). The predictions of LES 
are in good agreement with the value estimated from the experiment. These values 
are slightly lower than those obtained in the higher resolution LES with spectral 
eddy viscosity by Metais & Lesieur (1992). Results of running the simulation with 
no subgrid-scale model are also presented. 

The average energy spectrum E(k) is obtained in forced turbulence. Fig. 2 shows 
C K = e 2 / 3 k 5 / 3 E(k) plotted against the wavenumber k. According to Kolmogorov’s 
5/3 law (Kolmogorov, 1941), C K should be a constant in the inertial range. It is seen 
that the dynamic models with backscatter agree better with the 5/3 law than the 
purely dissipative models. Our best estimates for the Kolmogorov constant based on 
DLM(k) and DLM(S) are C K rs 1.8 and C K Rs 1.6 respectively. The experimentally 
measured values of C K are in the range 1.3 - 2.1 (Chasnov, 1991), though 1.5 is the 
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Figure 3. The level of backscatter as measured by energy transfer in forced 

isotropic turbulence at steady state : DLM(k); : DLM(S) and : 

as measured by fraction of points that have negative eddy-viscosity in DLM(k). 


commonly accepted value (Saddoughi & Veeravalli, 1994). The spectra predicted 
by the DM and DLM(+) are almost identical to each other, and they seem to decay 
somewhat faster than the 5/3 law. As expected these models without backscatter 
seem to be too dissipative. 

Fig. 3 shows the level of backscatter measured in two different ways. The solid 
line is the amount of energy being transferred from the subgrid to the large scales 
as a fraction of the net transfer els measured by 


/ [ r u 5 u. + ^ 3x 
1 / TijS t jd 3 x 


for the DLM(k) 


for the DLM(S) 

£sgs 

The dotted line is simply the fraction of grid points at which the Smagorinsky co- 
efficient is positive (only for the DLM(k); in the DLM(S) C is constrained to be 
positive). Here, we notice a substantial difference between the two models account- 
ing for backscatter. The deterministic DLM(k) predicts a much smaller amounts of 
backscatter than the stochastic DLM(S). 
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3. Future plans 

The first test of the DLM(S) has shown that this model is able to reproduce most 
of the feature of isotropic turbulence. However, the way backscatter is accounted for 
does not seem to play a major role in the present simulations. The main difference 
in the predictions obtained by these models concerns the rate of backscatter which 
is significantly higher in the DLM(S) than in the DLM(k). However, it may be 
unsafe to conclude that one of the models performs better from this difference. 
Indeed, there is not a lot of measurement of backscatter rate, and in addition, this 
rate has been shown to be filter- dependent in DNS (Piomelli et al 1991). Other 
performances of these two model are similar. 

DLM(S) is cheaper to implement, but the DLM(k) provides more information 
since it predicts the subgrid-scale energy as well as the pressure. Complex flows 
could differentiate further between these models, and the next step in the valida- 
tion of the stochastic model will be to implement it for more complex geometries. 
Preliminary work in the channel flow (Cabot, 1994) has shown that a slightly mod- 
ified version of the DLM(S) from which the pressure totally disappears could be 
implemented more easily than the formulation presented here. The advantage of 
developing a model without explicit need of the pressure is obvious in complex 
geometries in which the pressure has to be obtained by a u Poisson solver”. We 
thus plan to further develop this modified version of the DLM(S). Further tests in 
channel flows and mixing layer should follow. 
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